Instability of Bose-Einstein condensates in tilted lattices with time-periodical 

modulation 
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We study tiie dynamical stability of Bose-Einstein condensates in an optical lattice with a time- 
periodic modulation potential and a constant acceleration force simultaneously. We derive the 
explicit expressions of quasienergies and obtain the stability diagrams in the parameter space of the 
interaction strength and the modulation amplitude. The ratio of the acceleration force to the modu- 
lation frequency characterizes two cases: integer and non-integer resonances. For integer resonances, 
the critical interaction strength gc shows an alternate behavior where the completely unstable re- 
gions correspond to the negative effective tunneling strength. Among non-integer resonances, we 
observe that gc peaks are centered around half-integer resonances for which the completely un- 
stable regions disappear, accompanied with a whole displacement of gc- Compared with integer 
and half-integer resonances, the crossovers between them show no explicit dependence of gc on the 
modulation amplitude. 

PACS numbers: 67.85.-d, 03.65.Xp, 03.75.Kk, 03.75.Lm 



I. INTRODUCTION 

Since Bose-Einstein condensation (BEC) of dilute 
gases of alkali atoms was realized in 1995 Ir^], much 
attention has been paid to the ultracold atomic systems 
in various configurations, such as optical lattices which 
are formed by counter-propagating laser beams. Due to 
the flexibility of the system, the famous Mott-insulator- 
Superfluid transition of BECs in optical lattices has been 
observed and studied 4, 5]. Dynamically, there are also 
many interesting phenomena, such as the resonant tun- 
neling Q. In the limit of vanishing particle interactions, 
e.g., in dilute bosonic gases, the system may exhibit 
Bloch oscillation similar to the behavior of sin- 

gle electron in the crystalline field of solids. Other than 
the energetic instability 0) induced by the non-zero tem- 
perature, the nonlinear scattering process may also cause 
the dynamical instability 10] damping the oscillation be- 
havior 

Dynamical stability, which describes whether the sys- 
tem remains stable or not during its time evolution, is our 
concern in this paper. We are interested in the dynami- 
cal response of BECs in modulated optical lattices. We 
notice that dynamical stability of BECs system in an op- 
tical lattice modulated by a time-periodic potential has 
been investigated in both experiment and theory [l^] . 
Also, dynamical stability of the condensates in tilted op- 
tical lattices considering the effect of a constant accelera- 
tion force has attracted many people's attention [T^.[l5l|. 
Recently BECs in tilted time-periodically modulated lat- 
tices formed by both a time-periodic modulation poten- 
tial and a constant acceleration force was realized experi- 
mentally [IBl • We focus on the dynamical stability of the 
above system to which few attention has been paid be- 
fore, and aim to find out the influence of the interparticle 
interactions. 

The structure of this paper is organized as follows. In 



Sec.|TT]we describe our modulated optical lattice system 
with modified Bose-Hubbard Hamiltonian [l3, [l8|. In 
mean-field approximation we derive the time-evolution 
equations and obtain the explicit solutions. In order 
to determine the dynamical stability of the system, we 
adopt an ansatz of the time-evolution operator [Tsl. fl^. 
Upon Floquet theorem, the quasienergy analysis, which 
is strongly related to the dynamical stability, is given 
in Sec. IIIII Based on the ratio of the acceleration force 
to the modulation frequency, two specific cases, integer 
and non-integer resonances, are considered. We also dis- 
cuss the main features of dynamical stability for both 
two cases in Sec. IIVI A summary with brief discussion is 
given in Sec. IVl 



II. MODEL AND METHOD 

We consider ultracold bosonic gases on an one- 
dimensional optical lattice modulated by both a time- 
periodic potential and a constant external acceleration 
force simultaneously. To model such a system, we need 
to include the external modulation terms in the typical 
Bose-Hubbard Hamiltonian (BHH), i.e., 

H = - J ^ (aL^n + H.C.) + ^ ^ nr,iin„i - 1) 



-[Kcos{ujt) + d]Y^ 'I 



(1) 



Here Um (o-L) annihilates (creates) a boson on the mth 
lattice site, and hm = al^dm is the particle number oper- 
ator correspondingly. J describes the hopping strength 
between adjacent sites which are indicated by the sub- 
script (m, n). U refers to the on-site interaction between 
atoms. K and uj denote the amplitude and frequency of 
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the time-periodic modulation potential, while d stands 
for the constant acceleration force. 

The dynamical properties of our system can be ob- 
tained by studying the Heisenberg equations of motion 
for a„i. Provided that the particle number of atoms on 
each site is large enough, we can safely adopt the mean- 
field approximation (MFA) to replace these operators 
with their expectation values, e.g., Um = {dm) I ^/Nl, 
where is the average number of atoms per site. Then 
the time-evolution equations for are expressed as 



.dc 



dt 



= - J(am+1 + Ckm-l) + g\am P "r: 

+ [K cos.{ijjt) + d\mam, 



(2) 



where the natural unit ?i = 1 is taken and g = N^U 
denotes the normalized atomic interaction. Eqs. ^ are 
also regarded as discretized Gross-Pitaevskii (GP) equa- 
tions. Additionally, the particle number conservation 
Em^mOm = N gives \^rn? = L, wheiB N and L 
are the total numbers of atoms and optical lattice sites, 
respectively, and their ratio is Nl, i.e., Nl = N/L. 

In the interaction- free case, i.e., g — 0, Eqs. © de- 
scribe the single-particle dynamics whose analytical so- 
lutions can be explicitly derived through the gauge trans- 
formation |15i] 



am{t) exp 



then Eqs. ([2|) become 



dt 



K 



sin(a;i) 



am{t), (3) 



.darn 



m+1 t ^ L*m — 1 



(4) 



where D{t) = dt + K / uj sm{ujt) . Because the transla- 
tion symmetry is retrieved back in the above equation, 
we can impose the spatial periodic boundary conditions, 
i.e., am{t) — aL+m{t). The corresponding semi-classical 
Hamiltonian now reads 



H{t) = -Jj2{ 



(5) 

The Bloch-wave representation is desirable in order to 
obtain the solutions for am, 



~ r -1/2 „ikmu 

am = ^ ' 2^e bk. 



(6) 



where k — 2Trn/L is the quasimomentum (— tt < k < 
tt), n = 0,±1,...,±{L - l)/2 for odd L, while n = 
0, ±1, . . . , ±i/2 for even L. Then the evolution equa- 
tions of bk take the following simple form 



.dbk 



-2 J cos 



k-dt 



K 



sin(a;i) 



(7) 



whose solutions can be explicitly expressed as 

bu{t) = bkiO) exp {i2J [cos{k)S{t) - sin(fc)C(t)]} . (8) 



Here the factor 6fc(0) is the integral constant, and S{t) 
and C{t) are defined as 



S{t) 



sm{nuj + d)t 

71/, 1 -\- n 



+ CXD 



cos{nuj + d)t - 



V w 



(9) 



(10) 



where J7n denotes the nth ordinary Bessel function. 

The analytical solutions am for g = Q can be di- 
rectly obtained through the inverse Fourier transforma- 
tion. They are the starting points to derive the trial 
solutions for 5 7^ case. Specifically, we are concerned 
about 6fc(0) = (5fc_p case, which can simplify the final ex- 
pressions of solutions of Eqs. ^ as 



Oim{t) = exp{iTO[p — dt 



K 



sin(cjt)]} X 



exp{i2 J[cos(p)S'(i) - sin(p)C(t)] - igt}.{ll) 

Note that the solutions are applicable for L — cx) in 
which the boundary conditions make no difference. 

To explore the dynamical stability of the system, we 
assume small fluctuations around the stationary solutions 
following the usual ansatz [H, *-e., amit) = (^mi^) + 
(5a,„(t), where a^{t) take the expressions as Eqs. (fTTj) . 
The fluctuations can be expressed as 



5am{t) - a^(t)[u(t)e^«" + v* {t)e-"''' 



(12) 



Here q is the momentum of the excitation relative to the 
condensate. Substituting the trial solutions into Eqs. 
we obtain the Bogoliubov-de Gennes (BdG) equations for 
the quasiparticle excitations u{t) and v{t). 



^^(''^^=Mig,t)(''^'^ 



'dt\v{t)J '\v{t)J' 
where the elements of the matrix Ai{q, t) are given by 



(13) 



sin(cLit) — dt \ -\- g 



A^ii(q, i) = 4Jsin (^1^ sh " ^ 
Mi2{q,t) =g = -M2i{q,t), 

M22{q,t) = -4Jsin(^|^ sin ^| - ^ sin(wt) + dt ] - g. 

Note that the matrix M.{q,t) is also time-periodic with 
the periodicity T being the lowest common multiple of 
27r/w and iir/d. 

It will be convenient to introduce the evolution oper- 
ator U{t) in order to characterize the evolution of u{t) 
and w(t), i.e., {u{t),v{t) f = U{t){u{0),v{0))^ . Thus, 
the dynamical behavior U is governed by 



l-U{t)=M{q,t)U{t). 



(14) 



Using the 2x2 unit matrix as the initial data, we nu- 
merically solve Eq. (|14l) over period T. According to the 
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Floquet theorem, the eigenvalues of U(T) correspond 
to the excitation quasienergies via — exp[— ie^T] 
{i = 1,2). The dynamical stability of the system is spec- 
ified by the fact that both the quasienergies £i have no 
imaginary components, i.e., |Ai| = 1. If so for all values 
of q, then the solution is stable. Otherwise, the system 
may collapse. The analysis of those features can help 
us to map out the stability diagrams of the system in 
the parameter space of g/uj and K/uj which are shown in 
SecIVl 



IV. DYNAMICAL STABILITY 

We are interested in the dynamical stability of the sys- 
tem around the ground state, i.e., p — 0. To determine 
the critical interaction gc beyond which the system be- 
comes unstable, we increase the value of g from zero grad- 
ually for each K/uj until the unstable region is reached. 
Then we are able to plot the boundary between stable 
and unstable regions in the parameter space oi g/uo and 
K/oj. 



III. QUASIENERGY 

As stated above, quasienergies are tightly related to 
the dynamical stability. Since the period of the time- 
periodic evolution matrix in Eqs. @ is T, the Floquet 
theorem enables us to write the solution in the form of 
exp[— iet]'0(t), where e corresponds to the quasienergy 
of the system and 4i{t) corresponds to the Floquet state 
which shares the same period T. In order to obtain the 
quasienergies, we pick out all terms which are not T- 
periodic from a„i(i) p^ . Note that g = nuj+g' where n is 
a positive integer to ensure < g' < lo. Considering the 
possible divergence in the denominators of summations 
in S{t) and C{t), we analyze, separately, the cases of 
I = d/uj being integer and non-integer. 

For integer I, in the limit of d — > luj, both exp[sin(— ^ti;+ 
d)t/{—lu + d)] t and exp[— are not T-periodic, 
then we have the quasienergies expressed as 

e{p) = -2Jcos(p)J_( f - ) + g'. (15) 



2 4 6 8 10 12 14 2 4 6 8 10 12 14 
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FIG. 1. Quasienergy spectrum of a nine-site system for (a) 
I = 1 and (b) I — 2 and the K/lj dependences of the critical 
interaction Qc for (c) I = 1 and (d) 1 = 2. The displacement 
of quasienergy spectrum is g = g' = 0.5 for both cases. The 
vertical dashed lines are guides to eyes. 



The situation of d = has already been studied in 
Ref. so here we discuss more general case d ^ 0. 
Throughout our paper, J is set to be the unit. For exam- 
ple, we show the quasienergy spectrum for d = w (/ = 1) 
and d = 2uj {I = 2) in Fig. [UJa) and[Ub), respectively. 
The spectrum is obtained by the fundamental matrix 
method, which requires the diagonalization of q.{0)^ d{T) 
where d{t) — {ai{t),a2{t), . . . , aL{t)}. As L — >■ oo, the 
results match the analytical expressions ([TC)) perfectly. 
Another numerical way to obtain the quasienergies is to 
evolve the time-evolution operator of d{t) directly and 
diagonalize it, which yields the same results. It is obvi- 
ously observed from Fig.[lja) and[TIb) that the spectrum 
is wholly displaced by the amount of g' , as stated above. 
And the collapse points of the spectrum correspond to the 
zeroes of J'-i{K/ui), which manifests the phenomenon of 
the "coherent destruction of tunneling" (CDT) effect. 

For non- integer I, the case is much simpler since all 
terms in am{t) have the periodicity T except exp[— jg'i]. 
One can conclude that the quasienergies of the system 
are e{p) — g' , independent of the value of Kjui. Our 
numerical calculations further confirm our analysis. 



Since Eq. (fT4| can be solved analytically when 3 = 0, 
it is easy to find that the excitation quasienergies are all 
real such that the system is always stable regardless of 
the value of Kjbj. As g increases, the unstable regions 
emerge. We observe different stability diagrams for inte- 
ger and non- integer I. 

For the case of integer I, we plot the stability dia- 
grams for / = 1 and / = 2 in Fig. [IJc) and [Ijd) , re- 
spectively. It is clear that the dependence of g^ on Kjuj 
shows two kinds of alternative behaviors between plateau 
and upward concave shapes. We notice that the sep- 
arations between the two behaviors appear at the col- 
lapse points of the quasienergy spectrum, also the zeros 
of J-i{K/uj). From Eq. (|15p we may define the effective 
hopping strength Jcs{K/uj) = JJ^i{K/oj), while the fiat 
critical interaction regions (^c 0) correspond to the 
negative Jcg. The instability can be understood as the 
effect of attractive interactions [l^ for positive J^s due 
to the symmetry of the system Hamiltonian (g^l . When 
JceiK/io) = 0, CDT happens, freezing the dynamics of 
the system, thus insuring the stability, g^ peaks around 
these critical points. 
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Comparing Fig. [TJc) and [Hd) we notice a seemingly 
natural, however, interesting phenomenon. When I is a 
positive odd integer, in the limit of infinite small, however 
nonzero K/uj — >■ 0, the system is always unstable for any 
infinite small value of interaction, i.e., gc — >■ 0. The spe- 
cific choices of / insure the emergence of the phenomenon 
which has not been observed in the case of ^ = . 
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FIG. 2. The K/uj dependence of (a) Jo{K/uj) + J-xiKjuj) 
and (b) the critical interaction for I = 0.5. The vertical 
dashed lines are guides to eyes. 



are constant independent of Kjuj, no correspondence be- 
tween the quasienergies and peaks is observed, as ex- 
pected. In order to determine the positions of g^ peaks 
approximately, we have to resort to the expression of 
Cimify. Since the summations appear in the exponent of 
Eq. (jlll) . as |n-|-^| increases, the contribution from the n- 
th term reduces rapidly, if the high modulation frequency 
limit is assumed, i.e., a; 3> 1/J }21..] . The dominant terms 
would be those with the smallest values of |n -I- ?|. As an 
example, for I = 0.5, the dominant contributions to am(t) 
would come from Jq{K/uj) +J-i{K/uj) which is plotted 
in Fig. [2Ia). The correspondence between the zeros of 
Bessel functions and gc peaks is observed, confirming our 
analysis. As stated before, the dynamics of the system is 
frozen at the zeros of Bessel functions, indicating the ap- 
pearance of stable regions, where gc peaks with greatest 
probability. We also notice the whole displacement of gc 
whose value depends on d. In Fig.[3]we plot the behaviors 
of dynamical stability around integer and half-integer I, 
compared with which the crossover between them shows 
no explicit dependence of gc on the values of K/oj. We 
may interpret the behavior as the counterbalance of dif- 
ferent Bessel functions. Also from Fig. |31 we can find the 
same dependence of gc on d as found in Ref. ^15i] without 
time-periodical modulation, i.e., K/w — 0. 

Additionally, for both integer and non-integer Z, we 
find that the bigger value of w favors the high modulation 
frequency limit, also enhances the dynamical stability of 
the system as shown in Fig. [TJc) , [IJd) and Fig. [5l^b). 



SUMMARY 




FIG. 3. (Color online) Dependence of the critical interaction 
Qc on K and d for —1.5 < I < 1.5 in the case of cj = 4. 

For the case of non-integer I, the analysis is totally 
different, of which we start with half-integer I case. As 
an example, we plot the stability diagram for I = 0.5 
in Fig. m^b). Since the quasienergies for non-integer I 



We considered the system of BECs in one-dimensional 
tilted time-periodically modulated lattices. The ratio I of 
the acceleration force d to the modulation frequency uj is 
defined to help dividing our analysis into integer and non- 
integer resonances. Due to the time-periodicity of the 
system Hamiltonian, we gave the explicit expressions of 
quasienergies for both cases referring to integer and non- 
integer resonances by making use of Floquet theorem. 
For integer Z, the quasienergy spectrum shows its corre- 
spondence to the Bessel function J-i{K/uj) whose zeros 
correspond to the collapse points of the spectrum, indi- 
cating the appearance of CDT. However, for non-integer 
I, no dependence of the quasienergies on K/u) is found. 

We investigated the dynamical stability of the system 
by characterizing it in the parameter space of the in- 
teraction strength and the modulation amplitude. For 
integer Z, we observed two alternate behaviors of the de- 
pendence of the critical interaction strength gc on K/u: 
plateau {gc 0) and upward concave. The complete 
instability in horizontally fiat regions is brought out by 
the negative effective tunneling strength. The separa- 
tions between the two kinds of regions correspond to the 
zeros of Bessel Function J^i{K/uj). We also noticed an 
interesting phenomenon when I is a positive odd inte- 
ger. In such case, the system is always unstable in the 
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limit of K/uj — )■ corresponding to the infinite small, 
however, nonzero modulation potential. Among non- 
integer resonances, we observed that peaks are cen- 
tered around half-integer I for which the completely un- 
stable regions disappear. Also the whole displacement of 
gc{K/uj) emerges whose value is dependent on d. The po- 
sitions oi gc peaks in the direction of K/uj are determined 
approximately in the high modulation frequency limit. 
Compared with integer and half-integer resonances, the 
crossovers between them show no explicit dependence of 



f/c on the modulation amplitude. Since dynamical sta- 
bility is experimentally observable, these features are ex- 
pected to be confirmed in experiment. 
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